I. Getting raw data model-ready

Install Packages and set working directory

Read in Benthic Demography data and format

### "TCRMP_Master_Benthic_Cover_Feb2022.xlsx" was downloaded from the TCRMP data repository: https://sites.google.com/site/usvitcrmp/available-data?authuser=0
benthicFull <- read_excel("data/TCRMP_Master_Benthic_Cover_Feb2022.xlsx")%>%  
  rename(Site = Location, # rename Location as Site
         Date = FilmDate) %>% # This is date surveyed
  mutate(Date = as.Date(Date),  # format as date
         Year = year(Date), # make year column from survey date
         Month = month(Date), # make month column from survey date
         SampleID = paste(Site,Year,Month, sep="-")) %>% # SampleID column for unique survey points (all transects at site have same SampleID)
  relocate(c(SampleID, 
             Year, 
             Month), 
           .before= SampleYear) %>% #move important columns to front
  relocate(NoPts, .after = Date) %>% 
  select(-c(AnalysisBy, 
            AnalysisDate, 
            Period, 
            PC, 
            SampleYear, 
            SampleMonth)) #remove unused columns

Read in benthic codes data that matches TCRMP codes with functional groups

benthCodes<-read_xlsx("data/TcrmpBenthicCodes.xlsx") %>%
  filter(Code %in% names(benthicFull) )  # Subset to include only codes in benthicFull

Summarize data by functional group

#This collapses codes into total number of points for each category
benthCats <- benthicFull%>%
  pivot_longer(!c(SampleID, Site, Year, Transect, Month, Date, TWS, NoPts), names_to = "Code", values_to = "count") %>% #convert to long format for substitutions
  left_join(benthCodes, by="Code")%>% #join with codes reference table
  group_by(SampleID, Site, Year, Transect,Month, Date, NoPts, TWS, Category)%>% 
  summarise(count = sum(count)) %>%
  pivot_wider(id_cols = c(SampleID, Site, Year, Month, Date,Transect,TWS, NoPts), names_from = "Category", values_from = "count") %>% # pivot wider for categories as columns
  mutate(tottws = NoPts - TWS)%>%
  relocate(tottws, .before = "BLCov")%>%
  select(-c(TWS, UnkCov, OtherLiveCov, NonLiveCov, DisCov, BLCov)) #remove columns we won't use
## `summarise()` has grouped output by 'SampleID', 'Site', 'Year', 'Transect',
## 'Month', 'Date', 'NoPts', 'TWS'. You can override using the `.groups` argument.
## Adding missing grouping variables: `TWS`

Calculate % cover for each functional group

pctCov <- benthCats %>%
  mutate_at(vars(CalgCov:ZoanCov), 
            funs ((. / tottws))) %>%  
  rename_at(vars(CalgCov:ZoanCov), 
            list( ~paste("pctCov", gsub("_pctCov", "", .), sep = "_")))
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
## 
## # Simple named list: list(mean = mean, median = median)
## 
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
## 
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Import site metadata and add to pctCov

# Import site metadata
tcrmpMeta<-read_excel("data/TCRMP_SiteMetadata_SCTLDInf.xlsx")%>%
  rename(Site=Location)%>%
  mutate(InfectionDate = as.Date(InfectionDate, format =  "%m/%d/%Y"))

# Add site metadata to %cover dataframe
tcrmp_meta_pctcov <- pctCov%>%
  left_join(tcrmpMeta, by = "Site")%>%
  #relocate(c(Code:InfectionDate), .after = "Site")%>% 
  mutate(prePost = NA, #Make a pre-post column
  daysSinceInf=Date-InfectionDate) %>% #date minus infection date is number of days since infected (negative numbers tell you how soon before first sighting)
  relocate(Code:daysSinceInf, .after = "tottws")

Categorize each survey date as either pre or post infection based on infection dates

#populate the pre-post column
tcrmp_meta_pctcov$prePost[tcrmp_meta_pctcov$Date<=tcrmp_meta_pctcov$InfectionDate]<-"pre" #if date is before infection date, pre
tcrmp_meta_pctcov$prePost[tcrmp_meta_pctcov$Date>=tcrmp_meta_pctcov$InfectionDate]<-"post" #if date is after infection date, post

#-Make the days since variable numeric-#
tcrmp_meta_pctcov$daysSinceInf <- str_replace(tcrmp_meta_pctcov$daysSinceInf, " days", "")
tcrmp_meta_pctcov$daysSinceInf <- as.numeric(tcrmp_meta_pctcov$daysSinceInf)

#-Create a years since infection variable where all of the pre-SCTLD values are pooled-#
tcrmp_meta_pctcov$yearsSinceInf <- tcrmp_meta_pctcov$daysSinceInf/365
tcrmp_meta_pctcov$yearsSinceInf[which(tcrmp_meta_pctcov$yearsSinceInf < 0)] <- 0

tcrmp_meta_pctcov <- tcrmp_meta_pctcov %>%
  relocate(yearsSinceInf, 
           .after = daysSinceInf)

Save data to file

write.csv(tcrmp_meta_pctcov, "data/outputs/benthicTCRMP.csv")

Run brm of each transformed variable by “Years since infection”

First, set instructions for sampler

n_cores <- detectCores() # this will determine how many cores your computer has so that it can use all of its processing power
n_chains <- 2 
n_iter <- 8000
n_warmup <- 2000

#-Also add a function to make it easier when we're working with logit-transformed data-#
logit <- function(mu) {
  logit_mu <- log(mu/ (1 - mu))
  return(logit_mu)
}

#-Set inits-#
init_func <- function(chain_id = 1) {
  list ( Intercept  =  runif(1, .01, .1),
         b  =   runif(1, .001, .01),
         phi  =   runif(1, .001, .01),
         sd  =   runif(1, .001, .01),
         cor = runif(1, 0.001, .01),
         zi = runif(1, 0.001, .01))
}

# Create different inits for each chain

init_list <- vector(mode = "list", length = n_chains)

for(i in 1:n_chains)  
  {
  init_list[[i]] <- init_func(chain_id = i)
  }

CCA model

#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
  ggplot(aes(x = yearsSinceInf, y = pctCov_CalgCov)) +
  geom_point(aes(colour = Site)) +
  geom_smooth(method = "lm", se = F, aes(colour = Site)) +
  theme(legend.position = "none") + 
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

# subset just the Black Point site data because this will be reference level 
BlackPointData <- tcrmp_meta_pctcov %>%
  filter(Site=="Black Point")

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_CalgCov) #-we'll work with the zero-inflated beta

#define model
cca_mod <- bf(pctCov_CalgCov ~ 1  + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
              phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
              zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
              family =  zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))

##### Get priors
get_prior(cca_mod, data = tcrmp_meta_pctcov) 
##                 prior     class          coef group resp dpar nlpar lb ub
##                (flat)         b                                          
##                (flat)         b yearsSinceInf                            
##                lkj(1)       cor                                          
##                lkj(1)       cor                Site                      
##  student_t(3, 0, 2.5) Intercept                                          
##  student_t(3, 0, 2.5)        sd                                      0   
##  student_t(3, 0, 2.5)        sd                Site                  0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site                  0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site                  0   
##                (flat)         b                           phi            
##                (flat)         b yearsSinceInf             phi            
##  student_t(3, 0, 2.5) Intercept                           phi            
##  student_t(3, 0, 2.5)        sd                           phi        0   
##  student_t(3, 0, 2.5)        sd                Site       phi        0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site       phi        0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site       phi        0   
##                (flat)         b                            zi            
##                (flat)         b yearsSinceInf              zi            
##        logistic(0, 1) Intercept                            zi            
##  student_t(3, 0, 2.5)        sd                            zi        0   
##  student_t(3, 0, 2.5)        sd                Site        zi        0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site        zi        0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site        zi        0   
##        source
##       default
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_Calg <- BlackPointData %>%
              filter(pctCov_CalgCov> 0)

logit(mean(BPD_Calg$pctCov_CalgCov)) #-5.2
## [1] -5.213258
##### Set priors
cca_prior <- c(prior(normal(-5.2, 1), class = "Intercept"),
               prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
               prior(normal(0, 1), class = "b"),
               prior(normal(0, 1), class = "b", dpar = "zi"),
               prior(lkj(2), class = "cor"),
               prior(normal(0, 1), class = "Intercept", dpar = "phi"),
               prior(normal(0, 1), class = "b", dpar = "phi"),
               prior(normal(0, 1), class = "sd", dpar = "phi"),
               prior(exponential(1), class = "sd"),
               prior(exponential(1), class = "sd", dpar = "zi"))

##### Uncomment below to run model                         
# cca_brm <- brm(cca_mod,
#                data = tcrmp_meta_pctcov,
#                prior = cca_prior,
#                cores = n_cores,
#                chains = n_chains,
#                iter = n_iter,
#                warmup = n_warmup,
#                init = init_list)

##### Save model and/or read in saved model
#saveRDS(cca_brm, "data/outputs/cca_brm.RDS")
cca_brm <- readRDS("data/outputs/cca_brm.RDS")


##### Look at model       
cca_brm 
##  Family: zero_inflated_beta 
##   Links: mu = logit; phi = log; zi = logit 
## Formula: pctCov_CalgCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site) 
##          phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
##          zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
##    Data: tcrmp_meta_pctcov (Number of observations: 1330) 
##   Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
##          total post-warmup draws = 12000
## 
## Group-Level Effects: 
## ~Site (Number of levels: 34) 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                            0.75      0.10     0.58     0.98 1.00
## sd(yearsSinceInf)                        0.24      0.07     0.13     0.39 1.00
## sd(phi_Intercept)                        0.92      0.15     0.66     1.26 1.00
## sd(phi_yearsSinceInf)                    0.34      0.19     0.03     0.74 1.00
## sd(zi_Intercept)                         2.32      0.40     1.67     3.25 1.00
## sd(zi_yearsSinceInf)                     0.57      0.34     0.05     1.37 1.00
## cor(Intercept,yearsSinceInf)             0.37      0.23    -0.11     0.77 1.00
## cor(phi_Intercept,phi_yearsSinceInf)     0.22      0.34    -0.47     0.82 1.00
## cor(zi_Intercept,zi_yearsSinceInf)      -0.18      0.36    -0.82     0.54 1.00
##                                      Bulk_ESS Tail_ESS
## sd(Intercept)                            2770     5141
## sd(yearsSinceInf)                        5386     7242
## sd(phi_Intercept)                        4028     7082
## sd(phi_yearsSinceInf)                    2555     3744
## sd(zi_Intercept)                         3561     6221
## sd(zi_yearsSinceInf)                     2256     3399
## cor(Intercept,yearsSinceInf)             7530     8289
## cor(phi_Intercept,phi_yearsSinceInf)     9830     7529
## cor(zi_Intercept,zi_yearsSinceInf)       9834     7959
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept            -4.18      0.13    -4.44    -3.92 1.00     1655     2835
## phi_Intercept         4.92      0.18     4.56     5.27 1.00     3197     5215
## zi_Intercept         -2.01      0.44    -2.92    -1.17 1.00     2244     3352
## yearsSinceInf        -0.30      0.07    -0.43    -0.17 1.00     6367     7854
## phi_yearsSinceInf     0.38      0.14     0.11     0.65 1.00     8082     6806
## zi_yearsSinceInf      0.78      0.22     0.38     1.24 1.00     7904     7832
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#posterior predictive checks
pp_check(cca_brm, type = "hist", ndraws = 11) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(cca_brm, type = "boxplot", ndraws = 6)

pp_check(cca_brm, type = "dens_overlay", ndraws = 100) 

y <- posterior_predict(cca_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_CalgCov*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

Cyanobacteria model

#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
  ggplot(aes(x = yearsSinceInf, y = pctCov_CyanCov)) +
  geom_point(aes(colour = Site)) +
  geom_smooth(method = "lm", se = F, aes(colour = Site)) +
  theme(legend.position = "none") + 
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_CyanCov) #-we'll work with the zero-inflated beta

#define model
cyano_mod <- bf(pctCov_CyanCov ~ 1  + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
              phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
              zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
              family =  zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))

##### Get priors
get_prior(cyano_mod, data = tcrmp_meta_pctcov) 
##                 prior     class          coef group resp dpar nlpar lb ub
##                (flat)         b                                          
##                (flat)         b yearsSinceInf                            
##                lkj(1)       cor                                          
##                lkj(1)       cor                Site                      
##  student_t(3, 0, 2.5) Intercept                                          
##  student_t(3, 0, 2.5)        sd                                      0   
##  student_t(3, 0, 2.5)        sd                Site                  0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site                  0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site                  0   
##                (flat)         b                           phi            
##                (flat)         b yearsSinceInf             phi            
##  student_t(3, 0, 2.5) Intercept                           phi            
##  student_t(3, 0, 2.5)        sd                           phi        0   
##  student_t(3, 0, 2.5)        sd                Site       phi        0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site       phi        0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site       phi        0   
##                (flat)         b                            zi            
##                (flat)         b yearsSinceInf              zi            
##        logistic(0, 1) Intercept                            zi            
##  student_t(3, 0, 2.5)        sd                            zi        0   
##  student_t(3, 0, 2.5)        sd                Site        zi        0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site        zi        0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site        zi        0   
##        source
##       default
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_Cyano <- BlackPointData %>%
              filter(pctCov_CyanCov> 0)

logit(mean(BPD_Cyano$pctCov_CyanCov)) #-3.04
## [1] -3.043492
##### Set priors
cyano_prior <- c(prior(normal(-3, 1), class = "Intercept"),
               prior(normal(0, 1), class = "b"),
              prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
               prior(normal(0, 1), class = "b", dpar = "zi"),
               prior(exponential(1), class = "sd", dpar = "zi"),
               prior(lkj(2), class = "cor"),
               prior(normal(0, 1), class = "Intercept", dpar = "phi"),
               prior(normal(0, 1), class = "b", dpar = "phi"),
               prior(normal(0, 1), class = "sd", dpar = "phi"),
               prior(exponential(1), class = "sd"))

##### Uncomment below to run model              
# cyano_brm <- brm(cyano_mod,
#                data = tcrmp_meta_pctcov,
#                prior = cyano_prior,
#                cores = n_cores,
#                chains = n_chains,
#                iter = n_iter,
#                warmup = n_warmup,
#                init = init_list)

##### Save model and/or read in saved model
#saveRDS(cyano_brm, "data/outputs/cyanobacteria_brm.RDS")
cyano_brm <- readRDS("data/outputs/cyanobacteria_brm.RDS")

##### Look at model  
cyano_brm 
## Warning: There were 10 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: zero_inflated_beta 
##   Links: mu = logit; phi = log; zi = logit 
## Formula: pctCov_CyanCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site) 
##          phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
##          zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
##    Data: tcrmp_meta_pctcov (Number of observations: 1330) 
##   Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
##          total post-warmup draws = 12000
## 
## Group-Level Effects: 
## ~Site (Number of levels: 34) 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                            0.68      0.10     0.51     0.91 1.00
## sd(yearsSinceInf)                        0.71      0.13     0.48     0.99 1.00
## sd(phi_Intercept)                        0.88      0.13     0.66     1.17 1.00
## sd(phi_yearsSinceInf)                    0.47      0.27     0.03     1.02 1.00
## sd(zi_Intercept)                         0.93      0.16     0.65     1.27 1.00
## sd(zi_yearsSinceInf)                     1.63      0.49     0.84     2.72 1.00
## cor(Intercept,yearsSinceInf)             0.05      0.22    -0.37     0.47 1.00
## cor(phi_Intercept,phi_yearsSinceInf)    -0.43      0.34    -0.91     0.40 1.00
## cor(zi_Intercept,zi_yearsSinceInf)       0.07      0.22    -0.35     0.49 1.00
##                                      Bulk_ESS Tail_ESS
## sd(Intercept)                            2414     5154
## sd(yearsSinceInf)                        3001     6003
## sd(phi_Intercept)                        2954     4270
## sd(phi_yearsSinceInf)                    1070     2075
## sd(zi_Intercept)                         4594     6293
## sd(zi_yearsSinceInf)                     3536     5306
## cor(Intercept,yearsSinceInf)             2706     4261
## cor(phi_Intercept,phi_yearsSinceInf)     4957     5097
## cor(zi_Intercept,zi_yearsSinceInf)       3396     5732
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept            -3.57      0.12    -3.82    -3.32 1.00     1428     3003
## phi_Intercept         3.61      0.16     3.29     3.93 1.00     2005     3293
## zi_Intercept         -1.16      0.18    -1.53    -0.80 1.00     3122     6107
## yearsSinceInf         0.15      0.14    -0.14     0.43 1.00     2242     4272
## phi_yearsSinceInf     0.19      0.16    -0.10     0.55 1.00     2703     4353
## zi_yearsSinceInf     -0.86      0.38    -1.68    -0.17 1.00     4010     5783
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# posterior predictive checks
pp_check(cyano_brm, type = "hist", ndraws = 11) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(cyano_brm, type = "boxplot", ndraws = 6)

pp_check(cyano_brm, type = "dens_overlay", ndraws = 100)

y <- posterior_predict(cyano_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_CyanCov*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

Fire coral model

#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
  ggplot(aes(x = yearsSinceInf, y = pctCov_HydrozoanCoral)) +
  geom_point(aes(colour = Site)) +
  geom_smooth(method = "lm", se = F, aes(colour = Site)) +
  theme(legend.position = "none") + 
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_HydrozoanCoral) #-we'll work with zero-inflated beta

#define model
fireCoral_mod <- bf(pctCov_HydrozoanCoral ~ 1  + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
              phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
              zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
              family =  zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))

##### Get priors
get_prior(fireCoral_mod, data = tcrmp_meta_pctcov) 
##                 prior     class          coef group resp dpar nlpar lb ub
##                (flat)         b                                          
##                (flat)         b yearsSinceInf                            
##                lkj(1)       cor                                          
##                lkj(1)       cor                Site                      
##  student_t(3, 0, 2.5) Intercept                                          
##  student_t(3, 0, 2.5)        sd                                      0   
##  student_t(3, 0, 2.5)        sd                Site                  0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site                  0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site                  0   
##                (flat)         b                           phi            
##                (flat)         b yearsSinceInf             phi            
##  student_t(3, 0, 2.5) Intercept                           phi            
##  student_t(3, 0, 2.5)        sd                           phi        0   
##  student_t(3, 0, 2.5)        sd                Site       phi        0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site       phi        0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site       phi        0   
##                (flat)         b                            zi            
##                (flat)         b yearsSinceInf              zi            
##        logistic(0, 1) Intercept                            zi            
##  student_t(3, 0, 2.5)        sd                            zi        0   
##  student_t(3, 0, 2.5)        sd                Site        zi        0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site        zi        0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site        zi        0   
##        source
##       default
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_fireCoral <- BlackPointData %>%
              filter(pctCov_HydrozoanCoral> 0)

logit(mean(BPD_fireCoral$pctCov_HydrozoanCoral)) #-5.5
## [1] -5.450887
##### Set priors
fireCoral_prior <- c(prior(normal(-5.5, 1), class = "Intercept"),
               prior(normal(0, 1), class = "b"),
               prior(lkj(2), class = "cor"),
               prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
               prior(normal(0, 1), class = "b", dpar = "zi"),
               prior(exponential(1), class = "sd", dpar = "zi"),
               prior(normal(0, 1), class = "Intercept", dpar = "phi"),
               prior(normal(0, 1), class = "b", dpar = "phi"),
               prior(normal(0, 1), class = "sd", dpar = "phi"),
               prior(exponential(1), class = "sd"))

##### Uncomment to run model                         
# fireCoral_brm <- brm(fireCoral_mod,
#                data = tcrmp_meta_pctcov,
#                prior = fireCoral_prior,
#                cores = n_cores,
#                chains = n_chains,
#                iter = n_iter,
#                warmup = n_warmup,
#                init = init_list)

##### Uncomment below to save model and/or read in saved model
#saveRDS(fireCoral_brm, "data/outputs/fireCoral_brm.RDS")
fireCoral_brm <- readRDS("data/outputs/fireCoral_brm.RDS")

##### Look at model         
fireCoral_brm 
## Warning: There were 1 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: zero_inflated_beta 
##   Links: mu = logit; phi = log; zi = logit 
## Formula: pctCov_HydrozoanCoral ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site) 
##          phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
##          zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
##    Data: tcrmp_meta_pctcov (Number of observations: 1330) 
##   Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
##          total post-warmup draws = 12000
## 
## Group-Level Effects: 
## ~Site (Number of levels: 34) 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                            0.36      0.06     0.26     0.49 1.00
## sd(yearsSinceInf)                        0.07      0.05     0.00     0.19 1.00
## sd(phi_Intercept)                        0.79      0.14     0.56     1.11 1.00
## sd(phi_yearsSinceInf)                    0.17      0.14     0.01     0.53 1.00
## sd(zi_Intercept)                         1.55      0.24     1.15     2.08 1.00
## sd(zi_yearsSinceInf)                     0.32      0.23     0.01     0.86 1.00
## cor(Intercept,yearsSinceInf)             0.07      0.42    -0.73     0.81 1.00
## cor(phi_Intercept,phi_yearsSinceInf)    -0.05      0.40    -0.78     0.75 1.00
## cor(zi_Intercept,zi_yearsSinceInf)       0.05      0.41    -0.73     0.80 1.00
##                                      Bulk_ESS Tail_ESS
## sd(Intercept)                            4345     7094
## sd(yearsSinceInf)                        4756     5939
## sd(phi_Intercept)                        4246     6768
## sd(phi_yearsSinceInf)                    5087     6581
## sd(zi_Intercept)                         3683     6402
## sd(zi_yearsSinceInf)                     2769     5093
## cor(Intercept,yearsSinceInf)            13087     8056
## cor(phi_Intercept,phi_yearsSinceInf)    12654     8736
## cor(zi_Intercept,zi_yearsSinceInf)      13846     7685
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept            -5.21      0.08    -5.36    -5.06 1.00     3098     5274
## phi_Intercept         6.37      0.19     6.00     6.72 1.00     4121     5885
## zi_Intercept          0.98      0.29     0.41     1.55 1.00     2127     4044
## yearsSinceInf         0.08      0.05    -0.03     0.19 1.00     9591     9106
## phi_yearsSinceInf    -0.10      0.15    -0.41     0.20 1.00     9744     8595
## zi_yearsSinceInf     -0.22      0.15    -0.50     0.11 1.00     8079     5170
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# posterior predictive checks
pp_check(fireCoral_brm, type = "hist", ndraws = 11) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(fireCoral_brm, type = "boxplot", ndraws = 6)
## Notch went outside hinges
## ℹ Do you want `notch = FALSE`?
## Notch went outside hinges
## ℹ Do you want `notch = FALSE`?
## Notch went outside hinges
## ℹ Do you want `notch = FALSE`?
## Notch went outside hinges
## ℹ Do you want `notch = FALSE`?
## Notch went outside hinges
## ℹ Do you want `notch = FALSE`?
## Notch went outside hinges
## ℹ Do you want `notch = FALSE`?
## Notch went outside hinges
## ℹ Do you want `notch = FALSE`?

pp_check(fireCoral_brm, type = "dens_overlay", ndraws = 100)

y <- posterior_predict(fireCoral_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_HydrozoanCoral*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

Gorgonian model

#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
  ggplot(aes(x = yearsSinceInf, y = pctCov_GorgCov)) +
  geom_point(aes(colour = Site)) +
  geom_smooth(method = "lm", se = F, aes(colour = Site)) +
  theme(legend.position = "none") + 
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_GorgCov) #-we'll work with zero-inflated beta

#define model
gorgo_mod <- bf(pctCov_GorgCov ~ 1  + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
              phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
              zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
              family =  zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))

##### Get priors
get_prior(gorgo_mod, data = tcrmp_meta_pctcov) 
##                 prior     class          coef group resp dpar nlpar lb ub
##                (flat)         b                                          
##                (flat)         b yearsSinceInf                            
##                lkj(1)       cor                                          
##                lkj(1)       cor                Site                      
##  student_t(3, 0, 2.5) Intercept                                          
##  student_t(3, 0, 2.5)        sd                                      0   
##  student_t(3, 0, 2.5)        sd                Site                  0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site                  0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site                  0   
##                (flat)         b                           phi            
##                (flat)         b yearsSinceInf             phi            
##  student_t(3, 0, 2.5) Intercept                           phi            
##  student_t(3, 0, 2.5)        sd                           phi        0   
##  student_t(3, 0, 2.5)        sd                Site       phi        0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site       phi        0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site       phi        0   
##                (flat)         b                            zi            
##                (flat)         b yearsSinceInf              zi            
##        logistic(0, 1) Intercept                            zi            
##  student_t(3, 0, 2.5)        sd                            zi        0   
##  student_t(3, 0, 2.5)        sd                Site        zi        0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site        zi        0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site        zi        0   
##        source
##       default
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_gorgo <- BlackPointData %>%
              filter(pctCov_GorgCov> 0)

logit(mean(BPD_gorgo$pctCov_GorgCov)) #-3.6
## [1] -3.562578
##### Set priors
gorgo_prior <- c(prior(normal(-3.6, 1), class = "Intercept"),
               prior(normal(0, 1), class = "b"),
               prior(lkj(2), class = "cor"),
               prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
               prior(normal(0, 1), class = "b", dpar = "zi"),
               prior(exponential(1), class = "sd", dpar = "zi"),
               prior(normal(0, 1), class = "Intercept", dpar = "phi"),
               prior(normal(0, 1), class = "b", dpar = "phi"),
               prior(normal(0, 1), class = "sd", dpar = "phi"),
               prior(exponential(1), class = "sd"))

##### Uncomment to run model                
# gorgo_brm <- brm(gorgo_mod,
#                data = tcrmp_meta_pctcov,
#                prior = gorgo_prior,
#                cores = n_cores,
#                chains = n_chains,
#                iter = n_iter,
#                warmup = n_warmup,
#                init = init_list)

##### Save model and/or read in saved model
#saveRDS(gorgo_brm, "data/outputs/gorgonian_brm.RDS")
gorgo_brm<- readRDS("data/outputs/gorgonian_brm.RDS")

##### Look at model   
gorgo_brm 
##  Family: zero_inflated_beta 
##   Links: mu = logit; phi = log; zi = logit 
## Formula: pctCov_GorgCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site) 
##          phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
##          zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
##    Data: tcrmp_meta_pctcov (Number of observations: 1330) 
##   Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
##          total post-warmup draws = 12000
## 
## Group-Level Effects: 
## ~Site (Number of levels: 34) 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                            0.99      0.13     0.78     1.27 1.00
## sd(yearsSinceInf)                        0.09      0.05     0.01     0.21 1.00
## sd(phi_Intercept)                        0.99      0.15     0.73     1.32 1.00
## sd(phi_yearsSinceInf)                    0.14      0.10     0.01     0.39 1.00
## sd(zi_Intercept)                         2.69      0.47     1.92     3.75 1.00
## sd(zi_yearsSinceInf)                     0.39      0.25     0.02     0.94 1.00
## cor(Intercept,yearsSinceInf)            -0.23      0.35    -0.83     0.52 1.00
## cor(phi_Intercept,phi_yearsSinceInf)     0.08      0.41    -0.72     0.81 1.00
## cor(zi_Intercept,zi_yearsSinceInf)      -0.32      0.39    -0.90     0.58 1.00
##                                      Bulk_ESS Tail_ESS
## sd(Intercept)                            2856     5202
## sd(yearsSinceInf)                        4069     5645
## sd(phi_Intercept)                        4178     6217
## sd(phi_yearsSinceInf)                    5453     7333
## sd(zi_Intercept)                         3598     6354
## sd(zi_yearsSinceInf)                     4396     5902
## cor(Intercept,yearsSinceInf)            16965     8601
## cor(phi_Intercept,phi_yearsSinceInf)    19649     8670
## cor(zi_Intercept,zi_yearsSinceInf)      14832     9383
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept            -3.85      0.18    -4.21    -3.51 1.00     1305     2642
## phi_Intercept         4.66      0.18     4.30     5.03 1.00     2973     5193
## zi_Intercept         -2.52      0.52    -3.59    -1.54 1.00     2543     4046
## yearsSinceInf         0.01      0.04    -0.08     0.09 1.00    15958     9584
## phi_yearsSinceInf     0.03      0.10    -0.15     0.22 1.00    13837     9178
## zi_yearsSinceInf      0.08      0.22    -0.33     0.54 1.00    11903     8936
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#posterior predictive checks
pp_check(gorgo_brm, type = "hist", ndraws = 11) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(gorgo_brm, type = "boxplot", ndraws = 6)

pp_check(gorgo_brm, type = "dens_overlay", ndraws = 100)

y <- posterior_predict(gorgo_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_GorgCov*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

Macroalgae model

#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
  ggplot(aes(x = yearsSinceInf, y = pctCov_MacaCov)) +
  geom_point(aes(colour = Site)) +
  geom_smooth(method = "lm", se = F, aes(colour = Site)) +
  theme(legend.position = "none") + 
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_MacaCov) #-we'll work with beta

#define model
macro_mod <- bf(pctCov_MacaCov ~ 1  + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
              phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
                    family = Beta(link = "logit", link_phi = "log"))

##### Get priors
get_prior(macro_mod, data = tcrmp_meta_pctcov) 
##                 prior     class          coef group resp dpar nlpar lb ub
##                (flat)         b                                          
##                (flat)         b yearsSinceInf                            
##                lkj(1)       cor                                          
##                lkj(1)       cor                Site                      
##  student_t(3, 0, 2.5) Intercept                                          
##  student_t(3, 0, 2.5)        sd                                      0   
##  student_t(3, 0, 2.5)        sd                Site                  0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site                  0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site                  0   
##                (flat)         b                           phi            
##                (flat)         b yearsSinceInf             phi            
##  student_t(3, 0, 2.5) Intercept                           phi            
##  student_t(3, 0, 2.5)        sd                           phi        0   
##  student_t(3, 0, 2.5)        sd                Site       phi        0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site       phi        0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site       phi        0   
##        source
##       default
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_macro <- BlackPointData %>%
              filter(pctCov_MacaCov> 0)

logit(mean(BPD_macro$pctCov_MacaCov)) #-1
## [1] -0.9637347
##### Set priors
macro_prior <- c(prior(normal(-1, 1), class = "Intercept"),
               prior(normal(0, 1), class = "b"),
               prior(lkj(2), class = "cor"),
               prior(normal(0, 1), class = "Intercept", dpar = "phi"),
               prior(normal(0, 1), class = "b", dpar = "phi"),
               prior(normal(0, 1), class = "sd", dpar = "phi"),
               prior(exponential(1), class = "sd"))

##### Uncomment to run model
# macro_brm <- brm(macro_mod,
#                data = tcrmp_meta_pctcov,
#                prior = macro_prior,
#                cores = n_cores,
#                chains = n_chains,
#                iter = n_iter,
#                warmup = n_warmup,
#                init = init_list)

##### Save model and/or read in saved model
#saveRDS(macro_brm, "data/outputs/macroalgae_brm.RDS")
macro_brm<- readRDS("data/outputs/macroalgae_brm.RDS")

##### Look at model
macro_brm 
##  Family: beta 
##   Links: mu = logit; phi = log 
## Formula: pctCov_MacaCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site) 
##          phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
##    Data: tcrmp_meta_pctcov (Number of observations: 1330) 
##   Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
##          total post-warmup draws = 12000
## 
## Group-Level Effects: 
## ~Site (Number of levels: 34) 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                            0.89      0.11     0.69     1.14 1.00
## sd(yearsSinceInf)                        0.35      0.07     0.24     0.49 1.00
## sd(phi_Intercept)                        0.56      0.09     0.41     0.75 1.00
## sd(phi_yearsSinceInf)                    0.27      0.14     0.02     0.55 1.00
## cor(Intercept,yearsSinceInf)            -0.61      0.13    -0.82    -0.30 1.00
## cor(phi_Intercept,phi_yearsSinceInf)    -0.26      0.31    -0.77     0.47 1.00
##                                      Bulk_ESS Tail_ESS
## sd(Intercept)                            2026     3900
## sd(yearsSinceInf)                        3924     6727
## sd(phi_Intercept)                        4654     6889
## sd(phi_yearsSinceInf)                    2417     3204
## cor(Intercept,yearsSinceInf)             4609     6972
## cor(phi_Intercept,phi_yearsSinceInf)    10667     7094
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept            -0.71      0.15    -1.02    -0.41 1.00      797     1435
## phi_Intercept         3.11      0.11     2.90     3.32 1.00     3699     5961
## yearsSinceInf         0.12      0.07    -0.01     0.27 1.00     1711     3440
## phi_yearsSinceInf     0.48      0.10     0.29     0.69 1.00     9179     8812
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#posterior predictive checks
pp_check(macro_brm, type = "hist", ndraws = 11) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(macro_brm, type = "boxplot", ndraws = 6)

pp_check(macro_brm, type = "dens_overlay", ndraws = 100) #looks good

y <- posterior_predict(macro_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_MacaCov*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

SCTLD-resistant coral model

#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
  ggplot(aes(x = yearsSinceInf, y = pctCov_ResistantCoral)) +
  geom_point(aes(colour = Site)) +
  geom_smooth(method = "lm", se = F, aes(colour = Site)) +
  theme(legend.position = "none") + 
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_ResistantCoral) #-we'll work with zero-inflated beta

#define model
resistantCoral_mod <- bf(pctCov_ResistantCoral ~ 1  + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
              phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
              zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
              family =  zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))

##### Get priors
get_prior(resistantCoral_mod, data = tcrmp_meta_pctcov) 
##                 prior     class          coef group resp dpar nlpar lb ub
##                (flat)         b                                          
##                (flat)         b yearsSinceInf                            
##                lkj(1)       cor                                          
##                lkj(1)       cor                Site                      
##  student_t(3, 0, 2.5) Intercept                                          
##  student_t(3, 0, 2.5)        sd                                      0   
##  student_t(3, 0, 2.5)        sd                Site                  0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site                  0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site                  0   
##                (flat)         b                           phi            
##                (flat)         b yearsSinceInf             phi            
##  student_t(3, 0, 2.5) Intercept                           phi            
##  student_t(3, 0, 2.5)        sd                           phi        0   
##  student_t(3, 0, 2.5)        sd                Site       phi        0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site       phi        0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site       phi        0   
##                (flat)         b                            zi            
##                (flat)         b yearsSinceInf              zi            
##        logistic(0, 1) Intercept                            zi            
##  student_t(3, 0, 2.5)        sd                            zi        0   
##  student_t(3, 0, 2.5)        sd                Site        zi        0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site        zi        0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site        zi        0   
##        source
##       default
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_resistantCoral <- BlackPointData %>%
              filter(pctCov_ResistantCoral> 0)

logit(mean(BPD_resistantCoral$pctCov_ResistantCoral)) #-2.8
## [1] -2.807501
##### Set priors
resistantCoral_prior <- c(prior(normal(-2.8, 1), class = "Intercept"),
               prior(normal(0, 1), class = "b"),
               prior(lkj(2), class = "cor"),
               prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
               prior(normal(0, 1), class = "b", dpar = "zi"),
               prior(exponential(1), class = "sd", dpar = "zi"),
               prior(normal(0, 1), class = "Intercept", dpar = "phi"),
               prior(normal(0, 1), class = "b", dpar = "phi"),
               prior(normal(0, 1), class = "sd", dpar = "phi"),
               prior(exponential(1), class = "sd"))

##### Uncomment to run model                
# resistantCoral_brm <- brm(resistantCoral_mod,
#                data = tcrmp_meta_pctcov,
#                prior = resistantCoral_prior,
#                cores = n_cores,
#                chains = n_chains,
#                iter = n_iter,
#                warmup = n_warmup,
#                init = init_list)

##### Uncomment below to save model and/or read in saved model
#saveRDS(resistantCoral_brm, "data/outputs/resistantCoral_brm.RDS")
resistantCoral_brm <- readRDS("data/outputs/resistantCoral_brm.RDS")

##### Look at model       
resistantCoral_brm 
##  Family: zero_inflated_beta 
##   Links: mu = logit; phi = log; zi = logit 
## Formula: pctCov_ResistantCoral ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site) 
##          phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
##          zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
##    Data: tcrmp_meta_pctcov (Number of observations: 1330) 
##   Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
##          total post-warmup draws = 12000
## 
## Group-Level Effects: 
## ~Site (Number of levels: 34) 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                            0.71      0.09     0.56     0.91 1.00
## sd(yearsSinceInf)                        0.05      0.04     0.00     0.16 1.00
## sd(phi_Intercept)                        0.73      0.10     0.55     0.96 1.00
## sd(phi_yearsSinceInf)                    0.11      0.08     0.00     0.31 1.00
## sd(zi_Intercept)                         2.06      0.49     1.28     3.20 1.00
## sd(zi_yearsSinceInf)                     0.31      0.28     0.01     1.04 1.00
## cor(Intercept,yearsSinceInf)            -0.21      0.41    -0.84     0.69 1.00
## cor(phi_Intercept,phi_yearsSinceInf)    -0.16      0.43    -0.87     0.71 1.00
## cor(zi_Intercept,zi_yearsSinceInf)       0.07      0.44    -0.78     0.84 1.00
##                                      Bulk_ESS Tail_ESS
## sd(Intercept)                            1870     4118
## sd(yearsSinceInf)                        3458     5364
## sd(phi_Intercept)                        3240     5804
## sd(phi_yearsSinceInf)                    4628     5224
## sd(zi_Intercept)                         3423     5595
## sd(zi_yearsSinceInf)                     5332     5972
## cor(Intercept,yearsSinceInf)             8925     7147
## cor(phi_Intercept,phi_yearsSinceInf)    10718     8626
## cor(zi_Intercept,zi_yearsSinceInf)       9627     8870
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept            -3.59      0.12    -3.84    -3.35 1.00      903     1894
## phi_Intercept         4.69      0.14     4.42     4.95 1.00     2223     4119
## zi_Intercept         -4.48      0.55    -5.69    -3.54 1.00     3987     5574
## yearsSinceInf         0.02      0.03    -0.05     0.08 1.00     8267     7681
## phi_yearsSinceInf     0.05      0.08    -0.11     0.21 1.00    10306     9134
## zi_yearsSinceInf     -0.15      0.33    -0.86     0.45 1.00     8257     5633
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#posterior predictive checks
pp_check(resistantCoral_brm, type = "hist", ndraws = 11) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(resistantCoral_brm, type = "boxplot", ndraws = 6)

pp_check(resistantCoral_brm, type = "dens_overlay", ndraws = 100) 

y <- posterior_predict(resistantCoral_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_ResistantCoral*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

SCTLD-susceptible coral model

#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
  ggplot(aes(x = yearsSinceInf, y = pctCov_SuscCoral)) +
  geom_point(aes(colour = Site)) +
  geom_smooth(method = "lm", se = F, aes(colour = Site)) +
  theme(legend.position = "none") + 
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_SuscCoral) #-we'll work with zero-inflated beta

#define model
SuscCoral_mod <- bf(pctCov_SuscCoral ~ 1  + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
              phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
              zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
              family =  zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))

##### Get priors
get_prior(SuscCoral_mod, data = tcrmp_meta_pctcov) 
##                 prior     class          coef group resp dpar nlpar lb ub
##                (flat)         b                                          
##                (flat)         b yearsSinceInf                            
##                lkj(1)       cor                                          
##                lkj(1)       cor                Site                      
##  student_t(3, 0, 2.5) Intercept                                          
##  student_t(3, 0, 2.5)        sd                                      0   
##  student_t(3, 0, 2.5)        sd                Site                  0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site                  0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site                  0   
##                (flat)         b                           phi            
##                (flat)         b yearsSinceInf             phi            
##  student_t(3, 0, 2.5) Intercept                           phi            
##  student_t(3, 0, 2.5)        sd                           phi        0   
##  student_t(3, 0, 2.5)        sd                Site       phi        0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site       phi        0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site       phi        0   
##                (flat)         b                            zi            
##                (flat)         b yearsSinceInf              zi            
##        logistic(0, 1) Intercept                            zi            
##  student_t(3, 0, 2.5)        sd                            zi        0   
##  student_t(3, 0, 2.5)        sd                Site        zi        0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site        zi        0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site        zi        0   
##        source
##       default
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_SuscCoral <- BlackPointData %>%
              filter(pctCov_SuscCoral> 0)

logit(mean(BPD_SuscCoral$pctCov_SuscCoral)) #-2.4
## [1] -2.439212
##### Set priors
SuscCoral_prior <- c(prior(normal(-2.4, 1), class = "Intercept"),
               prior(normal(0, 1), class = "b"),
               prior(lkj(2), class = "cor"),
               prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
               prior(normal(0, 1), class = "b", dpar = "zi"),
               prior(exponential(1), class = "sd", dpar = "zi"),
               prior(normal(0, 1), class = "Intercept", dpar = "phi"),
               prior(normal(0, 1), class = "b", dpar = "phi"),
               prior(normal(0, 1), class = "sd", dpar = "phi"),
               prior(exponential(1), class = "sd"))

##### Uncomment to run model                       
# SuscCoral_brm <- brm(SuscCoral_mod,
#                data = tcrmp_meta_pctcov,
#                prior = SuscCoral_prior,
#                cores = n_cores,
#                chains = n_chains,
#                iter = n_iter,
#                warmup = n_warmup,
#                init = init_list)

##### Save model and/or read in saved model
#saveRDS(SuscCoral_brm, "data/outputs/susceptibleCoral_brm.RDS")
SuscCoral_brm <- readRDS("data/outputs/susceptibleCoral_brm.RDS")

##### Look at model
SuscCoral_brm 
##  Family: zero_inflated_beta 
##   Links: mu = logit; phi = log; zi = logit 
## Formula: pctCov_SuscCoral ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site) 
##          phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
##          zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
##    Data: tcrmp_meta_pctcov (Number of observations: 1330) 
##   Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
##          total post-warmup draws = 12000
## 
## Group-Level Effects: 
## ~Site (Number of levels: 34) 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                            1.12      0.15     0.88     1.46 1.00
## sd(yearsSinceInf)                        0.30      0.07     0.18     0.47 1.00
## sd(phi_Intercept)                        0.73      0.11     0.54     0.96 1.00
## sd(phi_yearsSinceInf)                    0.10      0.08     0.00     0.29 1.00
## sd(zi_Intercept)                         3.08      0.65     2.03     4.56 1.00
## sd(zi_yearsSinceInf)                     0.60      0.51     0.02     1.90 1.00
## cor(Intercept,yearsSinceInf)             0.05      0.23    -0.39     0.48 1.00
## cor(phi_Intercept,phi_yearsSinceInf)    -0.05      0.43    -0.81     0.77 1.00
## cor(zi_Intercept,zi_yearsSinceInf)       0.23      0.43    -0.69     0.88 1.00
##                                      Bulk_ESS Tail_ESS
## sd(Intercept)                            1651     2979
## sd(yearsSinceInf)                        3657     6269
## sd(phi_Intercept)                        3920     6756
## sd(phi_yearsSinceInf)                    5740     6482
## sd(zi_Intercept)                         4207     6759
## sd(zi_yearsSinceInf)                     3288     6503
## cor(Intercept,yearsSinceInf)             5678     8110
## cor(phi_Intercept,phi_yearsSinceInf)    15435     8051
## cor(zi_Intercept,zi_yearsSinceInf)      11588     8609
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept            -2.71      0.19    -3.09    -2.34 1.00      897     1550
## phi_Intercept         4.13      0.14     3.85     4.39 1.00     2833     5175
## zi_Intercept         -4.77      0.73    -6.35    -3.49 1.00     3206     5944
## yearsSinceInf        -0.59      0.07    -0.75    -0.45 1.00     4877     6782
## phi_yearsSinceInf     0.23      0.08     0.07     0.38 1.00    14689     9183
## zi_yearsSinceInf      0.71      0.39    -0.19     1.37 1.00     7712     6737
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# posterior predictive checks
pp_check(SuscCoral_brm, type = "hist", ndraws = 11) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(SuscCoral_brm, type = "boxplot", ndraws = 6)

pp_check(SuscCoral_brm, type = "dens_overlay", ndraws = 100)

y <- posterior_predict(SuscCoral_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_SuscCoral*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

Sponge model

#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
  ggplot(aes(x = yearsSinceInf, y = pctCov_SpoCov)) +
  geom_point(aes(colour = Site)) +
  geom_smooth(method = "lm", se = F, aes(colour = Site)) +
  theme(legend.position = "none") + 
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_SpoCov) #-we'll work with zero-inflated beta

#define model
sponge_mod <- bf(pctCov_SpoCov ~ 1  + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
              phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
              zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
              family =  zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))

##### Get priors
get_prior(sponge_mod, data = tcrmp_meta_pctcov) 
##                 prior     class          coef group resp dpar nlpar lb ub
##                (flat)         b                                          
##                (flat)         b yearsSinceInf                            
##                lkj(1)       cor                                          
##                lkj(1)       cor                Site                      
##  student_t(3, 0, 2.5) Intercept                                          
##  student_t(3, 0, 2.5)        sd                                      0   
##  student_t(3, 0, 2.5)        sd                Site                  0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site                  0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site                  0   
##                (flat)         b                           phi            
##                (flat)         b yearsSinceInf             phi            
##  student_t(3, 0, 2.5) Intercept                           phi            
##  student_t(3, 0, 2.5)        sd                           phi        0   
##  student_t(3, 0, 2.5)        sd                Site       phi        0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site       phi        0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site       phi        0   
##                (flat)         b                            zi            
##                (flat)         b yearsSinceInf              zi            
##        logistic(0, 1) Intercept                            zi            
##  student_t(3, 0, 2.5)        sd                            zi        0   
##  student_t(3, 0, 2.5)        sd                Site        zi        0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site        zi        0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site        zi        0   
##        source
##       default
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_sponge <- BlackPointData %>%
              filter(pctCov_SpoCov> 0)

logit(mean(BPD_sponge$pctCov_SpoCov)) #-2.4
## [1] -2.392563
##### Set priors
sponge_prior <- c(prior(normal(-2.4, 1), class = "Intercept"),
               prior(normal(0, 1), class = "b"),
               prior(lkj(2), class = "cor"),
               prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
               prior(normal(0, 1), class = "b", dpar = "zi"),
               prior(exponential(1), class = "sd", dpar = "zi"),
               prior(normal(0, 1), class = "Intercept", dpar = "phi"),
               prior(normal(0, 1), class = "b", dpar = "phi"),
               prior(normal(0, 1), class = "sd", dpar = "phi"),
               prior(exponential(1), class = "sd"))

##### Uncomment to run model                    
# sponge_brm <- brm(sponge_mod,
#                data = tcrmp_meta_pctcov,
#                prior = sponge_prior,
#                cores = n_cores,
#                chains = n_chains,
#                iter = n_iter,
#                warmup = n_warmup,
#                init = init_list)

##### Uncomment below to save model and/or read in saved model
#saveRDS(sponge_brm, "data/outputs/sponge_brm.RDS")
sponge_brm <- readRDS("data/outputs/sponge_brm.RDS")

##### Look at model               
sponge_brm 
##  Family: zero_inflated_beta 
##   Links: mu = logit; phi = log; zi = logit 
## Formula: pctCov_SpoCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site) 
##          phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
##          zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
##    Data: tcrmp_meta_pctcov (Number of observations: 1330) 
##   Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
##          total post-warmup draws = 12000
## 
## Group-Level Effects: 
## ~Site (Number of levels: 34) 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                            0.70      0.09     0.54     0.91 1.00
## sd(yearsSinceInf)                        0.18      0.05     0.10     0.29 1.00
## sd(phi_Intercept)                        0.41      0.08     0.26     0.59 1.00
## sd(phi_yearsSinceInf)                    0.28      0.11     0.07     0.52 1.00
## sd(zi_Intercept)                         2.63      0.61     1.64     4.01 1.00
## sd(zi_yearsSinceInf)                     0.52      0.43     0.02     1.61 1.00
## cor(Intercept,yearsSinceInf)             0.06      0.26    -0.44     0.56 1.00
## cor(phi_Intercept,phi_yearsSinceInf)     0.10      0.33    -0.51     0.75 1.00
## cor(zi_Intercept,zi_yearsSinceInf)       0.18      0.41    -0.67     0.86 1.00
##                                      Bulk_ESS Tail_ESS
## sd(Intercept)                            2750     4909
## sd(yearsSinceInf)                        4634     7231
## sd(phi_Intercept)                        5335     7016
## sd(phi_yearsSinceInf)                    3521     2731
## sd(zi_Intercept)                         4976     7898
## sd(zi_yearsSinceInf)                     4964     6747
## cor(Intercept,yearsSinceInf)             9675     8256
## cor(phi_Intercept,phi_yearsSinceInf)     8730     7172
## cor(zi_Intercept,zi_yearsSinceInf)      15921     9099
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept            -3.39      0.12    -3.63    -3.15 1.00     1221     2619
## phi_Intercept         4.72      0.09     4.55     4.90 1.00     6641     7590
## zi_Intercept         -5.10      0.70    -6.62    -3.90 1.00     6390     8149
## yearsSinceInf        -0.04      0.05    -0.15     0.06 1.00     6727     7309
## phi_yearsSinceInf    -0.02      0.11    -0.24     0.18 1.00    11167     8967
## zi_yearsSinceInf      0.43      0.41    -0.54     1.11 1.00    10812     7472
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#posterior predictive checks
pp_check(sponge_brm, type = "hist", ndraws = 11) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(sponge_brm, type = "boxplot", ndraws = 6)

pp_check(sponge_brm, type = "dens_overlay", ndraws = 100) 

y <- posterior_predict(sponge_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_SpoCov*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

Turf algae model

#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
  ggplot(aes(x = yearsSinceInf, y = pctCov_EACCov)) +
  geom_point(aes(colour = Site)) +
  geom_smooth(method = "lm", se = F, aes(colour = Site)) +
  theme(legend.position = "none") + 
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_EACCov) #-we'll work with beta

#define model
turf_mod <- bf(pctCov_EACCov ~ 1  + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
              phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
                    family = Beta(link = "logit", link_phi = "log"))

##### Get priors
get_prior(turf_mod, data = tcrmp_meta_pctcov) 
##                 prior     class          coef group resp dpar nlpar lb ub
##                (flat)         b                                          
##                (flat)         b yearsSinceInf                            
##                lkj(1)       cor                                          
##                lkj(1)       cor                Site                      
##  student_t(3, 0, 2.5) Intercept                                          
##  student_t(3, 0, 2.5)        sd                                      0   
##  student_t(3, 0, 2.5)        sd                Site                  0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site                  0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site                  0   
##                (flat)         b                           phi            
##                (flat)         b yearsSinceInf             phi            
##  student_t(3, 0, 2.5) Intercept                           phi            
##  student_t(3, 0, 2.5)        sd                           phi        0   
##  student_t(3, 0, 2.5)        sd                Site       phi        0   
##  student_t(3, 0, 2.5)        sd     Intercept  Site       phi        0   
##  student_t(3, 0, 2.5)        sd yearsSinceInf  Site       phi        0   
##        source
##       default
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_turf <- BlackPointData %>%
              filter(pctCov_EACCov> 0)

logit(mean(BPD_turf$pctCov_EACCov)) #-0.5
## [1] -0.5345477
##### Set priors
turf_prior <- c(prior(normal(-0.5, 1), class = "Intercept"),
               prior(normal(0, 1), class = "b"),
               prior(lkj(2), class = "cor"),
               prior(normal(0, 1), class = "Intercept", dpar = "phi"),
               prior(normal(0, 1), class = "b", dpar = "phi"),
               prior(normal(0, 1), class = "sd", dpar = "phi"),
               prior(exponential(1), class = "sd"))

##### Uncomment below to run model
# turf_brm <- brm(turf_mod,
#                data = tcrmp_meta_pctcov,
#                prior = turf_prior,
#                cores = n_cores,
#                chains = n_chains,
#                iter = n_iter,
#                warmup = n_warmup,
#                init = init_list)

##### Save model and/or read in saved model
#saveRDS(turf_brm, "data/outputs/turfAlgae_brm.RDS")
turf_brm<- readRDS("data/outputs/turfAlgae_brm.RDS")

##### Look at model
turf_brm 
##  Family: beta 
##   Links: mu = logit; phi = log 
## Formula: pctCov_EACCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site) 
##          phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
##    Data: tcrmp_meta_pctcov (Number of observations: 1330) 
##   Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
##          total post-warmup draws = 12000
## 
## Group-Level Effects: 
## ~Site (Number of levels: 34) 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                            0.65      0.08     0.51     0.84 1.00
## sd(yearsSinceInf)                        0.16      0.04     0.09     0.25 1.00
## sd(phi_Intercept)                        0.51      0.08     0.37     0.69 1.00
## sd(phi_yearsSinceInf)                    0.16      0.11     0.01     0.41 1.00
## cor(Intercept,yearsSinceInf)            -0.46      0.21    -0.81    -0.01 1.00
## cor(phi_Intercept,phi_yearsSinceInf)     0.02      0.38    -0.71     0.76 1.00
##                                      Bulk_ESS Tail_ESS
## sd(Intercept)                            2830     4796
## sd(yearsSinceInf)                        4595     7342
## sd(phi_Intercept)                        4712     6502
## sd(phi_yearsSinceInf)                    3214     6160
## cor(Intercept,yearsSinceInf)             9183     8281
## cor(phi_Intercept,phi_yearsSinceInf)    13776     7170
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept            -0.82      0.11    -1.04    -0.60 1.00     1511     2959
## phi_Intercept         3.29      0.10     3.09     3.49 1.00     5726     7532
## yearsSinceInf         0.03      0.04    -0.06     0.11 1.00     5649     7719
## phi_yearsSinceInf     0.10      0.10    -0.09     0.31 1.00     8783     8613
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#posterior predictive checks
pp_check(turf_brm, type = "hist", ndraws = 11) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(turf_brm, type = "boxplot", ndraws = 6)

pp_check(turf_brm, type = "dens_overlay", ndraws = 100)

y <- posterior_predict(turf_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_EACCov*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

Zoanthid model

  • omitted after conversation w/ T. Smith because he said they are extremely rare at the depths surveyed
#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
  ggplot(aes(x = yearsSinceInf, y = pctCov_ZoanCov)) +
  geom_point(aes(colour = Site)) +
  geom_smooth(method = "lm", se = F, aes(colour = Site)) +
  theme(legend.position = "none") + 
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_ZoanCov) #-we'll work with zero-inflated beta

ggplot(tcrmp_meta_pctcov)+
  geom_histogram(aes(x = pctCov_ZoanCov, fill = Site))+
                   facet_wrap(~Year)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

setup for generating predictions

#-Since our only two predictors in the model are Site and yearsSinceInf, this will let us generate predictions for all sites over a given period-#
SCTLD_sites <- unique(tcrmp_meta_pctcov$Site)
t <- seq(0, 10, by = (1/52)) #since effect is in years, this will estimate change per week

SCTLD_pred <- expand.grid(SCTLD_sites, t)

colnames(SCTLD_pred) <- c("Site", "yearsSinceInf")

Get R2 value for each model

bayes_R2(cca_brm)#.50
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.4992397 0.02361663 0.4518594 0.5452427
bayes_R2(cyano_brm)#.27
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.2652672 0.02774675 0.2129457 0.3217806
bayes_R2(gorgo_brm)#.67
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.6655684 0.01614981 0.6322491 0.6947094
bayes_R2(fireCoral_brm)#.24
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.2448967 0.03232055 0.1854114 0.3130079
bayes_R2(macro_brm)#.74
##     Estimate   Est.Error     Q2.5     Q97.5
## R2 0.7357125 0.005802696 0.723629 0.7464166
bayes_R2(resistantCoral_brm)#.48
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.4765436 0.03358091 0.4115251 0.5415283
bayes_R2(SuscCoral_brm)#.80
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.8005005 0.007277015 0.7848155 0.8131432
bayes_R2(sponge_brm)#.62
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.6215961 0.01437935 0.5916809 0.6483466
bayes_R2(turf_brm)#.65
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.6465213 0.009391979 0.6274814 0.6641061

get posterior effect estimates for plotting

posterior_draws_cca <- cca_brm %>%
 posterior_samples()%>%
  select(b_yearsSinceInf)%>%
  rename(effect_estimate = b_yearsSinceInf)%>%
  mutate(effect = "cca")%>%
  relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
posterior_draws_susc <- SuscCoral_brm %>%
 posterior_samples()%>%
  select(b_yearsSinceInf)%>%
  rename(effect_estimate = b_yearsSinceInf)%>%
  mutate(effect = "SCTLD-susceptible coral")%>%
  relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
posterior_draws_sponge <- sponge_brm %>%
 posterior_samples()%>%
  select(b_yearsSinceInf)%>%
  rename(effect_estimate = b_yearsSinceInf)%>%
  mutate(effect = "sponge")%>%
  relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
posterior_draws_resistant <- resistantCoral_brm  %>%
 posterior_samples()%>%
  select(b_yearsSinceInf)%>%
  rename(effect_estimate = b_yearsSinceInf)%>%
  mutate(effect = "SCTLD-resistant coral")%>%
  relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
posterior_draws_macroalg <- macro_brm %>%
 posterior_samples()%>%
  select(b_yearsSinceInf)%>%
  rename(effect_estimate = b_yearsSinceInf)%>%
  mutate(effect = "macroalgae")%>%
  relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
posterior_draws_turfalg <- turf_brm %>%
 posterior_samples()%>%
  select(b_yearsSinceInf)%>%
  rename(effect_estimate = b_yearsSinceInf)%>%
  mutate(effect = "turf algae")%>%
  relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
posterior_draws_gorgonian <- gorgo_brm %>%
 posterior_samples()%>%
  select(b_yearsSinceInf)%>%
  rename(effect_estimate = b_yearsSinceInf)%>%
  mutate(effect = "gorgo")%>%
  relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
posterior_draws_cyanobacteria <- cyano_brm %>%
 posterior_samples()%>%
  select(b_yearsSinceInf)%>%
  rename(effect_estimate = b_yearsSinceInf)%>%
  mutate(effect = "cyanobacteria")%>%
  relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
posterior_draws_fireCoral <- fireCoral_brm %>%
 posterior_samples()%>%
  select(b_yearsSinceInf)%>%
  rename(effect_estimate = b_yearsSinceInf)%>%
  mutate(effect = "fireCoral")%>%
  relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
PD_all <- posterior_draws_cca %>%
  full_join(posterior_draws_cyanobacteria)%>%
  full_join(posterior_draws_gorgonian)%>%
  full_join(posterior_draws_macroalg)%>%
  full_join(posterior_draws_resistant)%>%
  full_join(posterior_draws_sponge)%>%
  full_join(posterior_draws_susc)%>%
  full_join(posterior_draws_turfalg)%>%
  full_join(posterior_draws_fireCoral)%>%
  group_by(effect)%>%
  mutate(meanEffect = mean(effect_estimate))
## Joining with `by = join_by(effect, effect_estimate)`
## Joining with `by = join_by(effect, effect_estimate)`
## Joining with `by = join_by(effect, effect_estimate)`
## Joining with `by = join_by(effect, effect_estimate)`
## Joining with `by = join_by(effect, effect_estimate)`
## Joining with `by = join_by(effect, effect_estimate)`
## Joining with `by = join_by(effect, effect_estimate)`
## Joining with `by = join_by(effect, effect_estimate)`

# Fig 2A - plot of standardized effect size estimates

(tcrmpPosteriorDraws <- PD_all %>%
  ggplot(aes(x = effect_estimate, y = reorder(effect, -meanEffect))) +
  stat_halfeye()+
  geom_vline(xintercept = 0, size = 1, linetype = "dashed") +
  theme_bw(base_size = 20) +
  ylab("") +
  coord_cartesian(xlim = c(-1, 1))+
  xlab("Posterior estimate"))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

ggsave(file="fig/2a_tcrmpPosteriorDraws.svg", plot=tcrmpPosteriorDraws, width=10, height=10)

Figure 2B - plot of susceptible and resistant coral predictions

# FIGURE 2B
sitesIslands <- tcrmp_meta_pctcov%>%
  ungroup() %>%
  select(Site, Island) %>%
  unique()
#########################################Resistant coral #########################################

#-Predict changes to resistant at all sites over ten years-#
resistantCoral_pred <- predict(resistantCoral_brm, newdata = SCTLD_pred) #-make predictions from our new df-#

resistantCoral_pred <- bind_cols(SCTLD_pred, resistantCoral_pred) %>%
  as_tibble()%>%
  mutate(variable = "resistantCoral")%>%
  full_join(sitesIslands, by = "Site")

resistantCoral_site<-resistantCoral_pred%>%
  select(-variable)%>%
  group_by(Island, yearsSinceInf)%>%
  summarise(Q2.5 = mean(Q2.5), 
         Q97.5 = mean(Q97.5), 
         Estimate = mean(Estimate))
## `summarise()` has grouped output by 'Island'. You can override using the
## `.groups` argument.
#-Plot predicted resistant cover with 95% confidence intervals and real data-#
resistant10yrfig <- resistantCoral_site %>%
  ggplot(aes(x = yearsSinceInf, fill = Island)) +
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), alpha = 0.25) +
  geom_line(aes(y = Estimate, colour = Island), size = 2) +
  geom_point(data = tcrmp_meta_pctcov, aes(x = yearsSinceInf, y = pctCov_ResistantCoral, colour = Island)) +
  theme_classic(base_size = 20) + 
  ylim(0,0.25)+
 # theme(legend.position = "none", strip.background = element_blank()) +
  xlab("Years since infection") +
  ylab("resistant coral percent cover")+
  scale_fill_manual(values=c("#a1dab4", "#41b6c4", "#253494"))+
  scale_color_manual(values=c("#a1dab4", "#41b6c4", "#253494"))
resistant10yrfig
## Warning: Removed 4 rows containing missing values (`geom_point()`).

ggsave(file="fig/2b_resistant_10YearFig_Island.svg", plot=resistant10yrfig, width=10, height=10)
## Warning: Removed 4 rows containing missing values (`geom_point()`).
### MAKE RELATIVE DF FOR RESISTANT CORAL
resistantCoral_rel <- resistantCoral_pred %>%
  group_by(Site)%>%
  mutate(t0_resistant = Estimate[yearsSinceInf==0])%>%
  ungroup()%>%
  mutate(resistant_Rel = Estimate/t0_resistant)

# 5 year predicted mean change
resistantCoral_rel%>%
  filter(yearsSinceInf == 5)%>%
  summarise(pctChange = mean(resistant_Rel))
## # A tibble: 1 × 1
##   pctChange
##       <dbl>
## 1      1.11
# 10 year predicted mean change
resistantCoral_rel%>%
  filter(yearsSinceInf == 10)%>%
  summarise(pctChange = mean(resistant_Rel))
## # A tibble: 1 × 1
##   pctChange
##       <dbl>
## 1      1.32
#########################################Susceptible coral #########################################

#-Predict changes to Susc at all sites over ten years-#
SuscCoral_pred <- predict(SuscCoral_brm, newdata = SCTLD_pred) #-make predictions from our new df-#

SuscCoral_pred <- bind_cols(SCTLD_pred, SuscCoral_pred) %>%
  as_tibble()

SuscCoral_pred<-SuscCoral_pred%>%
  mutate(variable = "suscCoral")%>%
  full_join(sitesIslands, by = "Site")

suscCoral_site<-SuscCoral_pred%>%
  select(-variable)%>%
  group_by(Island, yearsSinceInf)%>%
  summarise(Q2.5 = mean(Q2.5), 
         Q97.5 = mean(Q97.5), 
         Estimate = mean(Estimate))
## `summarise()` has grouped output by 'Island'. You can override using the
## `.groups` argument.
#-Plot predicted resistant cover with 95% confidence intervals and real data-#
(susc10yrfig <- suscCoral_site %>%
  ggplot(aes(x = yearsSinceInf, fill = Island)) +
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), alpha = 0.25) +
  geom_line(aes(y = Estimate, colour = Island), size = 2) +
  geom_point(data = tcrmp_meta_pctcov, aes(x = yearsSinceInf, y = pctCov_SuscCoral, colour = Island)) +
  theme_classic(base_size = 20) + 
  ylim(0,0.25)+
 # theme(legend.position = "none", strip.background = element_blank()) +
  xlab("Years since infection") +
  ylab("resistant coral percent cover")+
  scale_fill_manual(values=c("#a1dab4", "#41b6c4", "#253494"))+
  scale_color_manual(values=c("#a1dab4", "#41b6c4", "#253494")))
## Warning: Removed 79 rows containing missing values (`geom_point()`).

ggsave(file="fig/2b_susc_10YearFig_Island.svg", plot=susc10yrfig, width=10, height=10)
## Warning: Removed 79 rows containing missing values (`geom_point()`).
### MAKE RELATIVE DF FOR SUSC CORAL

suscCoral_rel <- SuscCoral_pred %>%
  group_by(Site)%>%
  mutate(t0_susc = Estimate[yearsSinceInf==0])%>%
  ungroup()%>%
  mutate(susc_Rel = Estimate/t0_susc)

# 5 year predicted mean change
suscCoral_rel%>%
  filter(yearsSinceInf == 5)%>%
  summarise(pctChange = mean(susc_Rel))
## # A tibble: 1 × 1
##   pctChange
##       <dbl>
## 1    0.0821
# 10 year predicted mean change
suscCoral_rel%>%
  filter(yearsSinceInf == 10)%>%
  summarise(pctChange = mean(susc_Rel))
## # A tibble: 1 × 1
##   pctChange
##       <dbl>
## 1    0.0253